Orbital ice: an exact Coulomb phase on the diamond lattice 
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We demonstrate the existence of orbital Coulomb phase as the exact ground state of p-orbital ex- 
change Hamiltonian on the diamond lattice. The Coulomb phase is an emergent state characterized 
by algebraic dipolar correlations and a gauge structure resulting from local constraints (ice rules) of 
the underlying lattice models. For most ice models on the pyrochlore lattice, these local constraints 
are a direct consequence of minimizing the energy of each individual tetrahedron. On the contrary, 
the orbital ice rules are emergent phenomena resulting from the quantum orbital dynamics. We 
show that the orbital ice model exhibits an emergent geometrical frustration by mapping the de- 
generate quantum orbital ground states to the spin-ice states obeying the 2-in-2-out constraints on 
the pyrochlore lattice. We also discuss possible realization of the orbital ice model in optical lattices 
with p-band fermionic cold atoms. 

PACS numbers: 03.75.Ss, 05.50.+q, 71.10.Fd, 73.43.Nq 



I. INTRODUCTION 

Common water ice, a strongly correlated proton sys- 
tem, is a canonical example of geometrical frustration [l| . 
The oxygen ions in ice form a periodic diamond lattice 
whereas the protons are disordered due to the frustrated 
arrangement of two inequivalent O-H bonds with differ- 
ent lengths. This in turn leads to a macroscopic degener- 
acy of possible ground states and a finite entropy density 
of ice as temperature tends toward zero. Despite being 
disordered, the positioning of protons dictated by the so- 
called ice rules exhibits a strong short-range correlation 
in which each oxygen ion has two-near and two-far pro- 
tons. The ice rules also forbid single proton hopping and 
only allow for ring-exchange-type motion, reminiscent of 
the physics of gauge theory. 

A magnetic analogue of ice was discovered in py- 
rochlore oxides Dy2Ti2C>7 and Ho2Ti20y more than a 
decade ago @. These so-called spin ice compounds are 
essentially pyrochlore Ising magnet in which magnetic 
moments residing on a network of corner-sharing tetra- 
hedra are forced by single-ion anisotropy to point along 
the local (111) axes. It is found that extensively degener- 
ate spin configurations obeying the so-called '2-in-2-out' 
rules have essentially the same energy over a wide range 
of temperatures. The measured residual entropy is well 
approximated by the Pauling entropy for water ice §■ 
These local constraints require that every tetrahedron 
has two spins pointing in and two pointing out, in appar- 
ent analogy with the ice rules. Reversing a single spin 
in the ice state creates one defect tetrahedron with 3-in- 
1-out spins and another one with l-in-3-out spins. As 
recently pointed out in Ref. Q, these defect tetrahedra 
behave exactly as a gas of magnetic monopoles interact- 
ing with each other via Coulomb's 1/r law. 

Artificial versions of spin ice have also been created us- 
ing lithographically fabricated arrays of nanoscale mag- 
nets [a, l|| . Other proposals of artificial ice systems in- 
clude charged colloidals in optical traps and supercon- 
ducting vortices in specially fabricated pinning centers 



0, H| ■ A valence bond liquid phase with an ice- like degen- 
eracy is also shown to be the ground state of a spin-1/2 
Klein model on the pyrochlore lattice @. 

In most of these ice systems, the fundamental de- 
grees of freedom are doublet variables defined on the py- 
rochlore lattice or its two-dimensional counterpart. Their 
Hamiltonians can often be cast into the form 



eH' 



(1) 



where J > is the energy scale of excitations and the 
sum is over all tetrahedra. JC$$ is a nonnegative-definite 
operator defined for a tetrahedron. The last term denotes 
perturbations of energy scale e. The ice rules correspond 
to the contraints: 



/Cb = 0, 



(2) 



for all tetrahedra. Take spin ice as an example, the spin 
configurations can be specified by a set of Ising variables 
{(Ji\ such that Si = UiS e^, where e, denotes the lo- 
cal easy axis. The constraint operator is then given by 

ICei oc (Qk) 2 > where = J2 m m am is tne effective 
magnetic charge of a tetrahedron. The six up-up-down- 
down Ising configurations selected by constraints ^ cor- 
respond to the 2-in-2-out rules. Another example is the 
spin-1/2 Klein model for which IC^ = Vs m =2 is the pro- 
jection operator onto the subspace of maximum total spin 
= 2 9] . For temperatures in the regime f«r« J, 
configurations satisfying the 'ice rules' ([2]) for all tetra- 
hedra, are effectively degenerate. 

The ice model ([I} hosts an emergent Coulomb phase 
in which the local constraints K,^ = translate to a 
divergence-free flux V • B = in the coarse-grained 
approximation. The effective theory for the Coulomb 
phase is equivalent to conventional magnetostatics [To| . 
It follows that both the 'magnetic' field B and spins 
Si in this disordered yet highly constrained phase ex- 
hibit a dipolar-like correlation function (B a (0)Bp(r)) tx 
(5 a p — 3r a rp) /r 3 at large distances. 
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In this paper, we present an ice model in which the 
basic degrees of freedom are triplet orbital variables de- 
fined on the diamond lattice. Our investigation is partly 
motivated by recent progress in orbital-related many- 
body phenomena in optical lattices. We show that the 
strong directional dependence of orbital exchange com- 
bined with the special geometry of diamond lattice gives 
rise to a huge degeneracy in the Gutzwiller-type ground 
states which are also exact eigenstates of the orbital ex- 
change Hamiltonian. We demonstrate the existence of 
an orbital Coulomb phase by mapping the degenerate or- 
bital manifold to spin-ice states on the medial pyrochlorc 
lattice. It is worth noting that the orbital 'ice rules' are 
not constraints imposed by the Hamiltonian. Instead, 
they are emergent properties characterizing the short- 
range orbital correlations. This is in stark contrast to 
pyrochlore ice models (flj in which the ice rules are ex- 
plicitly incorporated in the Hamiltonian as the minimum 
energy states of individual tetrahedron. 



II. ORBITAL EXCHANGE HAMILTONIAN 

The ability to precisely control the interaction strength 
of cold atoms in optical lattices provides clean realiza- 
tions of strongly correlated models without many unde- 
sirable complexities usually encountered in material sys- 
tems [HI, [12]. I n particular, since the cold-atom systems 
are free of Jahn- Teller distortions, they offer a new op- 
portunity to investigate the intrinsic exchange physics 
associated with the orbital degrees of freedom [13]. One 
of the most interesting directions is the novel frustration 
phenomenon originating from the anisotropic orbital in- 
teractions. 

The exchange physics of p-orbitals in two-dimensional 
optical lattices has been extensively discussed in 
Refs. [3 Hq - The intricate interplay between lattice 
geometry and anisotropic orbital exchange leads to dra- 
matically distinct ground states in different lattices. The 
orbital exchange on a square lattice is dominated by an 
antiferromagnetic Ising-like Hamiltonian, which gives rise 
to a Neel-type orbital order. For triangular, honeycomb, 
and kagome lattices, the orbital interaction is described 
by a novel quantum 120° model. Although long-range or- 
bital orders occur in the cases of triangular and kagome 
lattices, orbital interactions are frustrated on the bipar- 
tite honeycomb lattice and a huge degeneracy remains 
in the classical ground state. These highly degener- 
ate ground states can be mapped to fully packed non- 
intersecting loops on the honeycomb lattice. Quantum 
fluctuations, on the other hand, select a six-site plaque- 
tte ground state through order from disorder mechanism. 
The 120° model also describes the effective orbital inter- 
action in transition metal oxides including honeycomb, 
cubic, and pyrochlore lattices [l6l - [l8j . 

Here we consider a p-band Hubbard model with spin- 
less fermions on three-dimensional optical lattices. We 
assume that each optical site is approximated by an 



isotropic harmonic potential. For two particles per site, 
one of them fills the inert s-orbital while the other one 
occupies one of the three p-orbitals. The kinetic terms 
of p-band fermions include a longitudinal tii and a trans- 
verse t±_ hopping, corresponding to a and 7r-bondings, 
respectively. Typically, £y ^> t± [19[ and we shall neglect 
the transverse hopping as a zeroth-order approximation. 
The fermions interact with each other through an on-site 
repulsion: H int = UJ2i, a ^/3 n ^ n ii3, where n ia = p\ a p la 
is the fermion number operator. The leading contribu- 
tion to U comes from the p-wave scattering for spinless 
fermions. The strong correlation regime U 3> tu can be 
potentially realized with the aid of the recently proposed 
stable optical p-wave Feshbach resonance [20] , which has 
the advantage of suppressing the high rate of three-body 
recombination. It should be noted that recent experi- 
mental results have shown some limitations of the opti- 
cal s-wave Feshbach scheme (2l| . Further experimental 
investigations are needed in order to verify the feasibility 
of increasing p-wave interaction through optical Feshbach 
resonance. 

With charge fluctuations suppressed in the Mott- 
insulating limit, there still remains a triplet orbital de- 
grees of freedom at each site. Exchange interactions 
between these localized orbital variables originate from 
the second-order virtual hopping of the fermions. Since 
we assume a dominating f ii , for a bond parallel to 
n = (n x ,n y ,n z ), longitudinal hopping is possible only 
when one of the particles occupies the orbital |n) = 



n x \p x ) 



\p y ) + n z \p z ), while the other one is in an 



orthogonal state. The energy gain of such an antiferro- 
orbital alignment is described by the Hamiltonian 



Hr,v — — J 



(3) 



(ij) 



Here J = t?,/U sets the exchange energy scale, / is the 

identity operator, and P nij = |ny)(fiy| is the projec- 
tion operator of the active orbital on a nearest-neighbor 
bond (ij). Obviously, the nature of the orbital exchange 
physics depends critically on the lattice geometry. 

III. CUBIC OPTICAL LATTICE 

As a warm-up, we first consider the case of cubic lat- 
tice [Fig. [Ha)]. Using a basis spanned by \p x ), \Py), and 
\p z ) states, the orbital projectors along the x, y, and z 
bonds can be expressed in terms of Gell-mann matrices 
A^ = diag(l,-l,0) and A^ = diag(l, 1, -2)/V3. By 
grouping them into a doublet operator r = (t x ,t v ) — 
(y / 3/2)(A (3) , A (8) ), the three orbital projectors are 

P a = (/ + 2r-e a )/3, {a = x,y,z), (4) 

with e x/y = (±^, |) and e z = (0,-1) [Fig. [Ha)]. The 
expectation value of the doublet vector (r) represents the 
disparities of on-site orbital occupation numbers. The 
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(a) 



(b) 



FIG. 1: (a) Cubic optical lattice, (b) Domain of doublet 
vector (t) for the cubic lattice. The three corners correspond 
to p x , p y , and p z orbitals, respectively. 



domain of (r) is an equilateral triangle [Fig.QJb)], whose 
three corners, (t) = e x , e y and e 2 , correspond to states 
with pure p x , p y , and p z orbitals, respectively. Substitut- 
ing the projectors P a into Eq. ^ , we obtain an effective 
Hamiltonian: 



-^cubic 



8 J 
9 



E 

a=x,y,z (ij)\\a 



( x * • g a) fa • e a ) 



(•5) 



up to an irrelevant constant cq = —AN J/ 3. Although 
Eq. (J5J) has the same form as the well-known 120° model, 
it is actually a classical Hamiltonian since the three or- 
bital projectors P a commute with each other. As a result, 
the eigenstates of -ffcubic are simultaneous eigenstates of 
the orbital occupation operators P" whose eigenvalues 
are or 1. Since each site has exactly one fermion, 
P x + P y + P z = 1, the orbital state at a given site can 
be specified by one of the three corners in the triangular 
domain of (t). Eq. ([5]) can then be viewed as a 3-state 
Potts model with anisotropic interactions. Take an x- 
bond for example, there are 3 different orbital configura- 
tions: (p x ,p x ), (p y / z ,Py/z), and (p x ,p y / z ) whose energies 
are 8J/9, 2J/9 and — 4J/9, respectively. 

To investigate the orbital correlations in the ground 
state, we performed classical Monte Carlo simulations 
with periodic boundary conditions on systems up to 
N = 24 3 sites. Figs. [5] (a) and (b) show the average 
bond energy e, specific heat c, and entropy density s 
as functions of temperature T. The bond energy ap- 
proaches eo = —2 J/9 asT^O, implying that 2/3 of the 
bonds with an energy e = — 4J/9 are in the antiferro- 
orbital ground states, while the remaining 1/3 are frus- 
trated with an energy of 2J/9. The macroscopic degen- 
eracy of the ground states is evidenced by a residual en- 
tropy density sq ~ 0.599 fcs obtained by integrating the 
specific- heat curve [Fig. HJb)]. The orbital correlation 
C T (r) = far) • t(0)) decays rather fast and is negligible 
beyond r ~ 5, indicating a disordered orbital liquid. At 
large separations, the correlation function decays expo- 
nentially as shown in Fig. [5] (c). 
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FIG. 2: Monte Carlo simulations of the classical Hamilto- 
nian ©. (a) and (b) show the temperature dependence of 
average bond energy e = {-ff C ubic)/3-/V, specific heat c and en- 
tropy density s, respectively. The dashed line in (b) indicates 
the entropy density In 3 at the high-temperature para-orbital 
phase. The orbital correlation function C T (r) = (r(r) ■ r(0)} 
is shown in (c) as a function of separation r. (d) shows 
lnCV(Z//2) as a function of linear system size L. 



The large residual entropy sq ps 0.599 also implies 
that the ground state is susceptible to nominally small 
perturbations present in the system. Indeed, as recently 
reported in Ref. 



25], 



inclusion of orbital interactions 
which break time-reversal symmetry induces long-range 
orbital ordering. As a final remark, it is worth noting 
that Eq. (J5J) is related to but quite different from the 
120° model with classical 0(2) spins, in which orbital- 
ordering is shown to be induced via order-from-disorder 
mechanism on the cubic lattice [22j-|24|. 



IV. DIAMOND OPTICAL LATTICE 

We now turn to orbital exchange on the oblique dia- 
mond lattice [Fig.[3ja)]. There are four distinct types of 
nearest-neighbor bonds pointing along directions fio = 
[111], n! = [111], n 2 = [111], and n 3 = [111]. Experi- 
mentally, a diamond optical lattice can be generated by 
the interference of four laser beams with a suitable ar- 
rangement of light polarizations 26] : 



V(r) 



x 



3 

E 

m— 1 



cos (K m • r) — cos (Kq • r) 



Here K m = (7r/2a)n m is the laser wave vector, and a 
is the nearest-neighbour bond length. To obtain the or- 



4 




(a) (b) 



FIG. 3: (a) Diamond optical lattice, (b) Domain of pseu- 
dovector (/x) for the diamond lattice. 

bital projectors on the nearest-neighbor bonds, we intro- 
duce a pseudovector fi = (fj, x , (j, y , p, z ) = (A^ 6 \ \( 4 \ A^) 
whose components are given by the three real-valued off- 
diagonal Gell-mann matrices. The operators /i a have the 
following nonzero elements: (p y \^ x \p z ) — (p z \fJ. y \p x ) — 
(j> x \n z \p y ) = 1. The orbital projectors along the four 
different bonds are 

P m = (/ + V3/x-n m )/3, (m = 0,l,2,3). (6) 

Substituting the above expression into Eq. ([3]) yields an 
effective Hamiltonian: 

2j 3 

-^diamond = — ^ ^ (/ij • n m ) (fij ■ n m ) . (7) 

m=0 (ij)\\m 

Since the three matrices \i a do not commute with each 
other, Eq. defines a quantum 'tetrahedral' Hamilto- 
nian for pseudovectors fii on the diamond lattice. The 
exchange interaction Q is geometrically frustrated. To 
see this, consider a bond (ij) along [111] direction. Its 
energy is minimized by orbital states \ipi) — \p x + p y + 
Pz)/VS and \ipj) — \p x ~ Py)/V%- The corresponding ex- 
pectation values of the pseudovector are (m) = 2h /y/3 
and (fij) = — z, respectively. However, such an antiferro- 
orbital alignment can not be achieved simultaneously on 
the other three (111) bonds attached to site i. 

In order to understand the ground-state structure, 
we first minimize the Hamiltonian using the Gutzwiller 
ansatz: 

i*)=ni^>=n>>^>- ^ 

i i 

The Gutzwiller wavefunction is a direct product of single- 
site orbitals, The orbital wavefunction at a given site is 
parameterized by two angles 9 and </>: 

\i/j) = sin 9 cos (f>\p x ) + sin 9 sin <j)\p y ) + cos 9\p z ) . 

The expectation value of the pseudovector is 

(/i) = (sin 2$ sin 0, sin 29 cos </>, sin 2 9 sin 2(f)) . (9) 



Fig. IS^b) shows the domain of (fi) which has a tetrahe- 
dral symmetry. We employ the Monte Carlo simulations 
to minimize the resulting mean- field energy E {(Hi)} = 
( v I'|i/diamond| x I'), which is a function of the pseudovectors. 
Specifically, small changes of 9i and <f>i are generated ran- 
domly and Eq. ^ is used to compute the change in (fii) 
and the corresponding AE. These updates are then ac- 
cepted according to detailed balancing. The Monte Carlo 
minimization yields many degenerate Gutzwiller ground 
states. We find that the pseudovectors in the ground 
states point along one of the six cubic directions, i.e., 
(fii) = ±x, ±y, or, ±z for all sites [Fig. 0], reminiscent 
of the six- vertex model. The corresponding orbital wave- 
functions are |±x) = \p y ±p z ) / y/2, and so on. The energy 
of each bond is exactly e — — 2J/9 in the ground state. 

Remarkably, the Gutzwiller ground states are also ex- 
act eigenstates of the Hamiltonian ([7]). To see this, we 
define an Ising variable for each of the nearest-neighbor 
bonds m attached to site i: 

aY 1 = V3{m)-n m = ±l, (m = 0,l,2,3). (10) 

They satisfy the orbital 'ice rules': 

<T?0? = -1 (11) 

for all nearest neighbors (ij) in the ground state. Now 
consider a given site i, if the Ising variable er™ = -1 on 
rn-th bond, is an eigenstate of the operator fii ■ n m 
with eigenvalue —l/y/3. On the other hand, for bonds 
with cr,- n = +1, an extra term is generated when acted 
by the same operator. Specifically, let \ipi) = |+x). The 
Ising variable is positive on [111] and [111] bonds; we 
have 

(m ■ n m ) |V0 = ±^V3\px) + 

with ± sign corresponding to m — and 1, respectively. 
Applying the combined bond operator on the Gutzwiller 
wavefunction yields 

(jii ■ n m ) (Mj ' n m ) |*) = TV2/3 \ Px ) t ® |f - 1/3 [*), 

where = Ylk^i \*Pk)- Note that the nearest-neighbor 
site j = j (to) depends on the bond index to. The two 
extra terms with opposite signs cancel each other when 
summed over m = and 1. The Gutzwiller state [VP) is 
thus an eigenstate of the sum of the two bond operators 
with positive cr|". Similar results hold for — \ ± y) 
or | ± z). Since each site has two bonds with a" 1 = +1 
attached to it, the extra terms cancel out when summed 
over all bonds. Consequently, the Gutzwiller state \^>) 
is an exact eigenstate of the full Hamiltonian. We also 
performed exact diagonalization of Eq. (|7|) on a finite 
system of 8 sites. With periodic boundary conditions, 
we find a huge degeneracy of the ground states which are 
indeed described by the Gutzwiller product. 

We now employ the fact that pyrochlore is the medial 
lattice of diamond to examine the degeneracy and struc- 
ture of the quantum ground states. As shown in Fig. 21 
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FIG. 4: A configuration of the pseudovectors on the diamond 
lattice and its mapping to the spin-ice state on the medial 
pyrochlore lattice. The pseudovector only assumes six dif- 
ferent values (fj,i) = ±x, ±y, and ±z in the ground states, 
corresponding to (p y ±p z ), (p z ±Px), and (p x ±p v ) orbitals, 
respectively. These six orbital configurations are mapped to 
the six 2-in-2-out ice states on a tetrahedron [Eq. (|12l) ]. 




2 4 6 8 10 12 14 0.5 1 1.5 2 2.5 

r InL 



FIG. 5: (a) Orbital correlation function C M (r) = {n{r) ■ fi(0)) 
as a function of distance r in the quantum ground state of 
Hamiltonian (b) shows lnC M (L/2) as a function of InL, 
here C M (L/2) is the correlation function between sites sep- 
arated by half the linear size L along a (110) chain of the 
lattice. The linear dependence in the log-log plot indicates a 
power-law decay: C M (L/2) ~ L~ 3 . 



V. SUMMARY AND DISCUSSION 



a pyrochlore magnet can be constructed by placing spins 
at the bond midpoints of a diamond lattice. This con- 
struction allows us to map the pseudovector field (fii) to 
a spin ice state on the pyrochlore lattice. Specifically, we 
label spins on a pyrochlore lattice by bond index (ij) of 
the diamond lattice and use the Ising variables (ITU)) to 
define its direction: 

S m =+^11™ = -of n m . (12) 

Here n m is a unit vector pointing from sites i to j. Note 
that the diamond-lattice sites are located at centers of 
tctrahcdra in the pyrochlore lattice, the above mapping 
shows that the six distinct values of pseudovectors in the 
ground state, i.e. (fii) = ±x, ±y, and ±z, correspond 
to the six different 2-in-2-out ice states on a tetrahedron 
as demonstrated in Fig. [4] The ground-state degeneracy 
of the diamond orbital model can thus be calculated us- 
ing the so-called Pauling estimate which gives a residual 
entropy per site so « £;s ln3/2 « 0.405 ks- 

The above mapping also makes it possible to com- 
pute orbital correlation functions by performing classical 
Monte Carlo simulations on pyrochlore spin ice. Since 
single-spin flip violates the ice rules, here we use the 
non-local loop moves to navigate the manifold of spin- 
ice ground states [27l.[28|: the results are shown in Fig. [5] 
The correlation function C^(r) = (fJ,(r) ■ fj,(0)) decays 
rather rapidly with the separation of spins. It is inter- 
esting to note that the pseudovector is related to the 
divergence- free flux via B(ri) ~ ±(/L»i), where ± sign 
refers to the two sublattices of the diamond lattice. As 
discussed in the introduction, the magnetic field B, hence 
the pseudovectors, display a dipolar-like correlation func- 
tion at long distances, as confirmed by our Monte Carlo 
simulations [Fig.[5fb)]. 



To summarize, we have investigated the orbital ex- 
change physics of p-band spinless fermions on both cu- 
bic and diamond lattices. In both cases we have found 
a macroscopic ground state degeneracy. The frustrated 
orbital interaction on the cubic lattice is governed by a 
classical three-state anisotropic Potts model. The ground 
state retains a finite entropy density Sq ~ 0.599fcg per 
site. Orbital correlation function decays exponentially at 
large distances. We have also derived a novel quantum 
'tetrahedral' model describing orbital interactions on the 
diamond lattice. We have obtained exact quantum many- 
body ground states which are extensively degenerate with 
a residual entropy density so ss ks In 3/2 0.405 ks- By 
mapping the degenerate quantum ground states to spin- 
ice states on a pyrochlore lattice, we have shown that 
the fermionic p-band Mott insulators on a diamond lat- 
tice can be viewed as an orbital analog of the frustrated 
ice phase I c of water. 

The huge degeneracy of orbital ice also helps circum- 
vent the entropy obstacle in its experimental realiza- 
tion. As noted in Ref. |29|, a major challenge in cre- 
ating strongly correlated phases in cold- atom systems is 
reaching the low level of entropies in such states. In this 
respect, the macroscopic residual entropy of the orbital 
ice renders the Coulomb phase much easier to realize in 
cold-atom optical lattices. 

It is worth noting that the orbital ice model presented 
in this paper is different in nature from most conventional 
ice systems. First, the fundamental degrees of freedom 
of orbital ice are orbital triplets defined on the diamond 
lattice, whereas those of the conventional ice models are 
Ising-like variables on pyrochlore. Second, the pyrochlore 
ice models with the ice rules explicitly incorporated into 
the Hamiltonian are essentially classical systems. On 
the other hand, the orbital ice is an intrinsic quantum 
model. The orbital 'ice rules' are emergent phenomena 
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resulting from the orbital exchange dynamics. This is a 
rare example of emergent geometrical frustration in three 
dimensions. As usually happens in highly frustrated sys- 
tems, the huge orbital degeneracy renders the ice phase 
susceptible to nominally small perturbations. Various in- 
teresting phases could emerge from the orbital Coulomb 



phase. Finally, it is also of great interest to examine the 
elementary excitations of the orbital ice model. 
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